Load packages

# Efficiently load packages 
pacman::p_load(tmap, terra, tmaptools, phyloseq, tidyverse, fun.gus, patchwork, sf, units, install = FALSE)

source("code/R/plotting_aesthetics.R")

load("data/07_biodiversity_exports/full_diversity_physeq.RData")

load("data/08_compositional_exports/full_abs_physeq.RData")

Load data

First, we need to connect our Rep_IDs with spatial information from our station log.

sample_sf <- full_diversity_physeq %>%
  sample_data() %>%
  data.frame() %>%
  st_as_sf(coords = c("Longitude","Latitude"), crs = "EPSG:4326")

# We use this object for surface interpolation
full_station_list_bio <- sample_sf %>%
  group_by(Station_ID) %>%
  slice_sample(n = 1) %>%
  ungroup() %>%
  select(Station_ID) %>%
  st_transform(crs = "EPSG:3174")

sample_sf_abs <- full_abs_physeq %>%
  sample_data() %>%
  data.frame() %>%
  st_as_sf(coords = c("Longitude","Latitude"), crs = "EPSG:4326")

I also need to read in a Lake Ontario outline I found from the USGS, here.

raw_ont <- read_sf("data/10_spatial_imports/hydro_p_LakeOntario/hydro_p_LakeOntario.shp") %>%
  st_transform(crs = "EPSG:4326")

ont_outline <- raw_ont %>%
  st_union()

raw_eri <- read_sf("data/10_spatial_imports/hydro_p_LakeErie/hydro_p_LakeErie.shp")

all_outlines <- list(raw_eri, raw_ont)

Reading in bathymetry

Downloaded my bathymetry from here

bathy <- terra::rast("data/10_spatial_imports/ontario_bath.tiff")

bathy_cropped <- mask(-bathy, vect(ont_outline))

I decided to use EPSG:3174 - NAD83/Great Lakes Albers. This is partially because sometime this data is combined with other Great Lakes data, and I want to use a projection that’s fairly good for all five lakes in the long term.

proj_bathy <- project(bathy_cropped, "EPSG:3174", method = "bilinear")

proj_ont <- st_transform(ont_outline, crs = "EPSG:3174") %>% st_sf()

Working with Enviromental Data from Metadata Exports

Reading in ctd data and converting to spatial object

load("data/05_metadata_exports/all_may_stations.RData")
load("data/05_metadata_exports/all_sep_stations.RData")
load("data/05_metadata_exports/clean_epa_chem.RData")

may_all_sp <- may_all_stations %>%
  group_by(Station_ID ) %>%
  summarize(Latitude = mean(latitude),
            Longitude = mean(longitude)) %>% 
  st_as_sf(coords = c("Longitude", "Latitude"), crs = "EPSG:4326") %>%
  ungroup()

sep_all_sp <- sep_all_stations %>%
  group_by(Station_ID) %>%
  summarize(Latitude = mean(latitude),
            Longitude = mean(longitude)) %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = "EPSG:4326") %>%
  ungroup()

This next section is a little tricky. We want to show our transects as lines, and we’ll use this for our vertical cross sections of temperature (Figure S5). Everything below is me makes these nice transect lines and extracting bathymetry along those transects.

shoreline_spots_both <- sep_all_sp %>%
  filter(Station_ID%in% c(8,17,29, 35, 38, 43, 48, 717, 62, 66)) %>%
  mutate(Station_ID = Station_ID + .5) %>%
  st_nearest_points(ont_outline) %>%
  st_cast("POINT") %>%
  st_transform(crs = "EPSG:3174")
  
shoreline_spots <- st_sf(Station_ID = .5 + c(8,17,29, 35, 38, 43, 48, 62,66, 717),
        geometry = shoreline_spots_both[seq(2,20,2)])

may_lines <- may_all_sp %>%
  filter(!Station_ID %in% c(71,74)) %>%
  st_transform(crs = "EPSG:3174") %>%
  rbind(shoreline_spots) %>%
  mutate(Transect = case_when(Station_ID %in%c(8:19, 8.5, 17.5)~"T1",
                            Station_ID %in%c(27:35, 29.5, 35.5)~"T2",
                            Station_ID %in%c(38:45, 43.5, 38.5)~"T3",
                            Station_ID %in%c(48:58, 716, 717, 48.5, 717.5)~"T4",
                            Station_ID %in% c(62:67,715, 66.5, 62.5)~"T5")) %>% 
  mutate(lat = st_coordinates(.)[,2]) %>%
  arrange(Transect, lat) %>% 
  group_by(Transect) %>%
  summarize(Stations = n(),
            do_union = FALSE) %>%
  st_cast(to = "LINESTRING") %>%
  st_transform(crs = "EPSG:3174")


samples <- st_line_sample(may_lines, sample = seq(from = 0, to = 1, length.out = 100)) %>% st_sf() %>% st_cast("POINT")

lengths <- st_length(may_lines)

sep_lines <- sep_all_sp %>%
  filter(!Station_ID %in% c(71,74)) %>%
  st_transform(crs = "EPSG:3174") %>%
  rbind(shoreline_spots) %>%
  mutate(Transect = case_when(Station_ID %in%c(8:19, 8.5, 17.5)~"T1",
                            Station_ID %in%c(27:35, 29.5, 35.5)~"T2",
                            Station_ID %in%c(38:45, 43.5, 38.5)~"T3",
                            Station_ID %in%c(48:58, 716, 717, 48.5, 717.5)~"T4",
                            Station_ID %in% c(62:67,715, 66.5, 62.5)~"T5")) %>%
  mutate(lat = st_coordinates(.)[,2]) %>%
  arrange(lat) %>% 
  group_by(Transect) %>%
  summarize(Stations = n(),
            do_union = FALSE) %>%
  st_cast(to = "LINESTRING") %>%
  st_transform(crs = "EPSG:3174")

samples_sep <- st_line_sample(sep_lines, sample = seq(from = 0, to = 1, length.out = 100)) %>% st_sf() %>% st_cast("POINT")
lengths_sep <- st_length(sep_lines)

tm_shape(proj_ont) + 
  tm_borders() + 
tm_shape(sep_lines) +
  tm_lines(col = "Transect",
           col.scale = tm_scale_categorical(values = "set1"),
           col.legend = tm_legend(title = "",
                                  text.size = 12,
                                  position = tm_pos_out("right","center"),
                                  frame = FALSE))+
  tm_layout(frame = FALSE)

We also map the GLATOs buoy locations on top of these transect lines, for the inset map in Figure S6A.

buoys <- data.frame(buoy = c("LNR314","WLO003","PPW016", "LON024", "LOJ038", "LOH024", "Lake Erie at Buffalo"),
           lat = c(43.3408, 43.51568, 43.77682, 43.41419,43.67258,43.81907, 42.887067),
           lon = c(-79.08577,-79.44942, -77.17267, -77.7472,-76.43503,-77.72727, -78.908267)) %>%
  st_as_sf(coords = c("lon", "lat"), crs = "EPSG:4326") %>%
  st_transform(crs = "EPSG:3174")

st_bbox(buoys)
##      xmin      ymin      xmax      ymax 
## 1404091.9  717197.3 1644997.4  827839.1
uni_outlines <- map(all_outlines, st_union)

uni_proj <- map(uni_outlines, st_transform, crs = "EPSG:3174")

com_outlines <- st_combine(do.call("c", uni_proj))

tm_shape(com_outlines, bbox = st_bbox(buoys) * c(.96, .98, 1.02, 1.08)) + 
  tm_borders() + 
  tm_shape(sep_lines) + 
  tm_lines(col = "grey") + 
  tm_shape(buoys) + 
  tm_symbols(fill = "black", size = .5, shape = 21
            ) +
  tm_text(text = "buoy", size = 1, ymod = 1, xmod = 1, col.legend = tm_legend(show = FALSE)) + 
  tm_compass(position = c(0,.8), size = 2, text.size = 1) + 
  tm_scalebar(breaks = c(0,25, 50), position = c(.8,0), text.size = 1)

Here, we extract the bathymetry along our transects in May and September. We also transformed our CTD casts and sampling points to “Distance Along Transect”, rather than Lat/Long.

# Extract bathymetry
may_bath_values <- terra::extract(proj_bathy, samples, xy = TRUE) %>%
  mutate(Transect = rep(c("T1","T2","T3","T4","T5"), each = 100),
         index_along_transect = rep(seq(from = 0, to = 1, length.out = 100), 5),
         transect_length = rep(lengths, each = 100),
         dist_along_transect = as.numeric(index_along_transect * transect_length)) %>%
  dplyr::select(ontario_bath, Transect, dist = dist_along_transect) %>%
  mutate(ontario_bath = ifelse(is.na(ontario_bath), 0, ontario_bath))

may_station_dists <- may_all_sp %>%
  filter(!Station_ID %in% c(71,74)) %>%
  st_transform(crs = "EPSG:3174") %>%
  rbind(shoreline_spots) %>%
  mutate(Transect = case_when(Station_ID %in%c(8:19, 8.5, 17.5)~"T1",
                            Station_ID %in%c(27:35, 29.5, 35.5)~"T2",
                            Station_ID %in%c(38:45, 43.5, 38.5)~"T3",
                            Station_ID %in%c(48:58, 716, 717, 48.5, 717.5)~"T4",
                            Station_ID %in% c(62:67,715, 66.5, 62.5)~"T5")) %>% 
  mutate(x = st_coordinates(.)[,1],
         y = st_coordinates(.)[,2]) %>%
  group_by(Transect) %>%
  arrange(y, .by_group = TRUE) %>%
  mutate(lag_x = lag(x),
         lag_y = lag(y),
         dist = sqrt((lag_x - x)^2 + (lag_y - y)^2),
         dist = ifelse(is.na(dist),0,dist),
         cum_dist = cumsum(dist)) %>%
  dplyr::select(Station_ID, dist = cum_dist, Transect)

may_with_dist <- may_all_stations %>%
  left_join(may_station_dists)



sep_bath_values <- terra::extract(proj_bathy, samples_sep, xy = TRUE) %>%
  mutate(Transect = rep(c("T1","T2","T3","T4","T5"), each = 100),
         index_along_transect = rep(seq(from = 0, to = 1, length.out = 100), 5),
         transect_length = rep(lengths_sep, each = 100),
         dist_along_transect = as.numeric(index_along_transect * transect_length)) %>%
  dplyr::select(ontario_bath, Transect, dist = dist_along_transect)

sep_station_dists <- sep_all_sp %>%
  filter(!Station_ID %in% c(71,74)) %>%
  st_transform(crs = "EPSG:3174") %>%
  rbind(shoreline_spots) %>%
  mutate(Transect = case_when(Station_ID %in%c(8:19, 8.5, 17.5)~"T1",
                            Station_ID %in%c(27:35, 29.5, 35.5)~"T2",
                            Station_ID %in%c(38:45, 43.5, 38.5)~"T3",
                            Station_ID %in%c(48:58, 716, 717, 48.5, 717.5)~"T4",
                            Station_ID %in% c(62:67,715, 66.5, 62.5)~"T5")) %>%
  mutate(x = st_coordinates(.)[,1],
         y = st_coordinates(.)[,2]) %>%
  group_by(Transect) %>%
  arrange(y, .by_group = TRUE) %>%
  mutate(lag_x = lag(x),
         lag_y = lag(y),
         dist = sqrt((lag_x - x)^2 + (lag_y - y)^2),
         dist = ifelse(is.na(dist),0,dist),
         cum_dist = cumsum(dist)) %>%
  dplyr::select(Station_ID, dist = cum_dist, Transect)

sep_with_dist <- sep_all_stations %>%
  left_join(sep_station_dists)

sep_sam_points <- sample_sf_abs %>%
  as.data.frame() %>%
  filter(month == "September")%>%
  select(depth = Deployment_Depth_m, Transect = transect, Station_ID, Upwelling) %>%
  left_join(sep_station_dists)

may_sam_points <- sample_sf_abs %>%
  as.data.frame() %>%
  filter(month == "May")%>%
  select(depth = Deployment_Depth_m, Transect = transect, Station_ID, Upwelling) %>%
  left_join(may_station_dists)

Finally, we’ll inteprolat and plot those vertical transects, each as its own little panel, which we’ll combine together with patchwork.

trans <- c("T1","T2","T3","T4","T5")

ODV_colours <- c( "#d31f2a", "#ffc000", "#0db5e6", "#7139fe")

temp_limits <- range(c(may_with_dist$temperature, sep_with_dist$temperature)) * c(0.7, 1.2)

temp_limit_may <- range(may_with_dist$temperature) * c(0.7, 1.2)

depth_limit <- max(c(may_with_dist$depth, sep_with_dist$depth)) + 10

title_cols <- RColorBrewer::brewer.pal(5, "Set1")

may_plots <- map2(trans, title_cols, function(trans,title_col){
  t1 <- may_with_dist %>%
    dplyr::filter(Transect == trans) %>%
    mutate(dist = dist/1000)
  
  ctd_mba <- MBA::mba.surf(t1[c("dist", "depth", "temperature")], no.X = 300, no.Y = 300, extend = T, sp = TRUE, b.box = c(0,75,0,300), h = 4)
  
  inter_temp <- data.frame(ctd_mba$xyz.est@coords, ctd_mba$xyz.est@data) %>%
    dplyr::rename(dist = x, depth = y, temperature = z)

  bottoms <- may_bath_values %>%
    dplyr::filter(Transect == trans) %>%
    dplyr::select(depth = ontario_bath, dist) %>%
    rbind(data.frame(depth = c(depth_limit, depth_limit),
                      dist = c(max(.$dist), min(.$dist)))) %>%
    mutate(dist = dist / 1000)
  

  ggplot(inter_temp, aes(x = dist, y = depth, fill = temperature)) +
    geom_raster() +
    geom_contour(aes(z = temperature), breaks = 6, color = "white", linetype = 2) + 
    geom_line(data = t1, aes(group = Station_ID), alpha = 1)  +
    geom_polygon(data = bottoms, fill = "grey40") +
    geom_point(data = filter(may_sam_points, Transect == trans), aes(y = depth, x = dist/1000, shape = Upwelling), color = "black", fill = "white", size = 2.5) + 
    scale_shape_manual(values = upwelling_shapes, guide = "none") + 
    scale_y_reverse() +
    scale_fill_gradientn(colors = rev(ODV_colours), limits = temp_limit_may) + 
    coord_cartesian(xlim = range(bottoms$dist), 
                    ylim = c(depth_limit,0),
                    expand = FALSE) +
    labs(x = "Distance Along Transect (km)", y = "Depth (m)", fill = "Temperature (°C)", title = paste0("May Transect: ", trans)) +
    theme(plot.title = element_text(color = title_col))
})
 
sep_plots <- map2(trans, title_cols, function(trans, title_col){
  t1 <- sep_with_dist %>%
    dplyr::filter(Transect == trans) %>%
    mutate(dist = dist/1000)
  
  ctd_mba <- MBA::mba.surf(t1[c("dist", "depth", "temperature")], no.X = 300, no.Y = 300, extend = T, sp = TRUE, b.box = c(0,75,0,300), h = 4)
  
  inter_temp <- data.frame(ctd_mba$xyz.est@coords, ctd_mba$xyz.est@data) %>%
    dplyr::rename(dist = x, depth = y, temperature = z)
  
  bottoms <- sep_bath_values %>%
    dplyr::filter(Transect == trans) %>%
    dplyr::select(depth = ontario_bath, dist) %>%
    rbind(data.frame(depth = c(depth_limit, depth_limit),
                     dist = c(max(.$dist), min(.$dist)))) %>%
    mutate(dist = dist / 1000)
 
  ggplot(inter_temp, aes(x = dist, y = depth, fill = temperature)) +
    geom_raster() + 
    geom_contour(aes(z = temperature), breaks = 12, color = "white", linetype = 2) + 
    geom_line(data = t1, aes(group = Station_ID), alpha = 1) + 
    geom_polygon(data = bottoms, fill = "grey40") + 
    geom_point(data = filter(sep_sam_points, Transect == trans), aes(y = depth, x = dist/1000, shape = Upwelling), color = "black", fill = "white", size = 2.5) + 
    scale_shape_manual(values = upwelling_shapes, guide = "none") + 
    scale_y_reverse() +
    scale_fill_gradientn(colors = rev(ODV_colours), limits = temp_limits)+
    coord_cartesian(xlim = range(bottoms$dist), ylim = c(max(bottoms$depth), 0 ), expand = FALSE) + 
    labs(x = "Distance Along Transect (km)", y = "Depth (m)", fill = "Temperature (°C)", title = paste0("Sept Transect: ", trans)) +
    theme(plot.title = element_text(color = title_col))
})

all_plots <- wrap_plots(may_plots, guides = "collect", ncol = 5) / wrap_plots(sep_plots, ncol = 5, guides = "collect")


all_plots

Writing out data for interpolation in QGIS

In the next sections, we write out geopackages which we then visualize in QGIS.

Writing out data for nutrient interpolation

surface <- rbind(may_all_stations,sep_all_stations) %>%
  group_by(Station_ID, month) %>%
  mutate(max_depth = max(depth)) %>%
  dplyr::filter(abs(depth - 5) == min(abs(depth - 5))) %>%
  slice_min(order_by = depth) %>%
  ungroup() %>%
  st_as_sf(coords = c("longitude","latitude"), crs = "EPSG:4326") %>%
  st_transform(crs = "EPSG:3174") %>%
  left_join(
    filter(clean_epa_chem, Depth_Class == "E")
    ) %>%
  mutate(TN_TP = TN/TP * 2.21)

full_station_list <- surface %>%
  select(Station_ID) %>% 
  group_by(Station_ID) %>%
  slice_sample(n = 1)

# Plotting Nutrients
split_nutrients <- surface %>%
  select(Station_ID, month, NH4:TN_TP) %>%
  group_by(month) %>%
  group_split()

st_write(split_nutrients[[1]], "analysis/QGIS_Work/may_all_nutrients.gpkg", append = FALSE)
## Deleting layer `may_all_nutrients' using driver `GPKG'
## Writing layer `may_all_nutrients' to data source `analysis/QGIS_Work/may_all_nutrients.gpkg' using driver `GPKG'
## Writing 31 features with 17 fields and geometry type Point.
st_write(split_nutrients[[2]], "analysis/QGIS_Work/sept_all_nutrients.gpkg", append = FALSE)
## Deleting layer `sept_all_nutrients' using driver `GPKG'
## Writing layer `sept_all_nutrients' to data source `analysis/QGIS_Work/sept_all_nutrients.gpkg' using driver `GPKG'
## Writing 30 features with 17 fields and geometry type Point.

Writing out data for cell count interpolation

may_cells <- sample_sf %>%
  filter(month=="May", Depth_Class == "E") %>%
  select(Richness, avg_cells_per_ml) %>%
  st_transform(crs = "EPSG:3174")

 st_write(may_cells, 
          "analysis/QGIS_Work/May_Cells.gpkg",
          append = FALSE)
## Deleting layer `May_Cells' using driver `GPKG'
## Writing layer `May_Cells' to data source `analysis/QGIS_Work/May_Cells.gpkg' using driver `GPKG'
## Writing 15 features with 2 fields and geometry type Point.
sept_cells <- sample_sf %>%
  filter(month=="September", Depth_Class == "E") %>%
  select(Richness, avg_cells_per_ml) %>%
  st_transform(crs = "EPSG:3174")

st_write(sept_cells, 
          "analysis/QGIS_Work/September_Cells.gpkg", 
          append = FALSE)
## Deleting layer `September_Cells' using driver `GPKG'
## Writing layer `September_Cells' to data source `analysis/QGIS_Work/September_Cells.gpkg' using driver `GPKG'
## Writing 15 features with 2 fields and geometry type Point.
cell_count_sp <- sample_sf %>%
  st_transform(crs = "EPSG:3174")

Plotting Harmful Algae (Figure S13)

melted <- full_abs_physeq %>%
  psmelt()

dol_df <- melted %>%
  select(Rep_ID, Genus, Species, Abundance, Comp_Group_Hier) %>%
  filter(str_detect(Genus, "Dolichospermum")) %>% 
  group_by(Rep_ID, Genus, Comp_Group_Hier) %>%
  summarize(Dol_Abund = sum(Abundance)) %>% ungroup()

micro_df <- melted %>%
  select(Rep_ID, Genus, Species, Abundance, Comp_Group_Hier) %>%
  filter(str_detect(Genus, "Microcystis")) %>% 
  group_by(Rep_ID, Genus, Comp_Group_Hier) %>%
  summarize(Mic_Abund = sum(Abundance)) %>% ungroup()

dol_df %>% 
  ggplot(aes(x = Comp_Group_Hier, y = Dol_Abund, color = Comp_Group_Hier)) + 
  geom_boxplot(outliers = FALSE) + 
  ggbeeswarm::geom_beeswarm() + 
  ggpubr::stat_compare_means(comparison = list(c(1,2), c(1,3), c(2,3)),
                             label = "p.signif") + 
  scale_color_manual(values = comp_three_colors) + 
  scale_x_discrete(labels = c("Deep", "Shallow\nMay","Shallow\nSeptember")) +
    scale_y_continuous(expand = expansion(mult = c(0.1,.1)))+
  theme(legend.position = "none") + 
  labs(x = "", y = "Dolichospermum NIES41 sp.\n(cells/mL)")

micro_df %>% 
  ggplot(aes(x = Comp_Group_Hier, y = Mic_Abund, color = Comp_Group_Hier)) + 
  geom_boxplot(outliers = FALSE) + 
  ggbeeswarm::geom_beeswarm() + 
  ggpubr::stat_compare_means(comparison = list(c(1,2), c(1,3), c(2,3)),
                             label = "p.signif") + 
  scale_color_manual(values = comp_three_colors) + 
  scale_x_discrete(labels = c("Deep", "Shallow\nMay","Shallow\nSeptember")) +
  scale_y_continuous(expand = expansion(mult = c(0.1,.1)),
                     labels = scales::label_comma())+
  theme(legend.position = "none") + 
  labs(x = "", y = "Microcystis PCC-7914 sp.\n(cells/mL)")

Then we write out values for spatial interpolation in QGIS

harmful_sp <- dol_df %>%
  left_join(cell_count_sp) %>%
  filter(month == "September", Depth_Class == "E") %>%
  left_join(select(micro_df, Rep_ID, Mic_Abund)) %>%
  select(Rep_ID, Dol_Abund, Mic_Abund, geometry)

st_write(harmful_sp, "analysis/QGIS_Work/harmful_algae.gpkg", append = FALSE)
## Deleting layer `harmful_algae' using driver `GPKG'
## Writing layer `harmful_algae' to data source `analysis/QGIS_Work/harmful_algae.gpkg' using driver `GPKG'
## Writing 15 features with 3 fields and geometry type Point.

Distance Decays and Average Dissimilarities

This section uses distance matrices defined in 08_Compositional_Analysis to test relationships between distance and depth decay, as well as variance partitioning.

load("data/08_compositional_exports/absolute_wunifrac.RData")

abs_unifrac_dist_0.5 <- absolute_wunifrac$unifracs[, , "d_0.5"]

unweighted_unifrac_dist <- absolute_wunifrac$unifracs[, , "d_UW"]
diag(abs_unifrac_dist_0.5) <- NA # We don't want identical comparisons

mat_df <- as.data.frame(abs_unifrac_dist_0.5) # Convert to data.frame

# Pivot to long format where each row is a pairwise comparison
tidy_dist_wun <- mat_df %>%
    mutate(item1 = rownames(mat_df)) %>%
    pivot_longer(!item1, names_to = "item2", values_to = "dist") %>%
  filter(!is.na(dist))

# Repeat for unweighted unifrac

diag(unweighted_unifrac_dist) <- NA

mat_df_un <- as.data.frame(unweighted_unifrac_dist)

tidy_dist_un <- mat_df_un %>%
    mutate(item1 = rownames(mat_df_un)) %>%
    pivot_longer(!item1, names_to = "item2", values_to = "un_dist") %>%
  filter(!is.na(un_dist))

# Join our distances

tidy_dist <- inner_join(tidy_dist_wun, tidy_dist_un)

# Select metadata for comparisons

rep_to_comp <- sample_data(full_abs_physeq) %>%
    data.frame() %>%
    dplyr::select(Rep_ID,
           Comp_Group_Hier,
           Station_ID,
           Latitude,
           Longitude,
           month,
           Depth_Class)
  
# Find average dissimilarities for each sample within it's Comp Group

avg_dists <- tidy_dist %>%
    mutate(
      item1_group = rep_to_comp$Comp_Group_Hier[match(item1, rep_to_comp$Rep_ID)],
      item2_group = rep_to_comp$Comp_Group_Hier[match(item2, rep_to_comp$Rep_ID)]
    ) %>%
    dplyr::filter(item1_group == item2_group) %>%
    group_by(item1) %>%
    summarize(avg.dist = mean(dist, na.rm = TRUE),
              avg.un_dist =mean(un_dist, na.rm = TRUE)) %>%
    rename(Rep_ID = item1)

We write out these values, for possible visualization in QGIS (I don’t think that made it into the paper tho)

split_dists <- 
  avg_dists %>% 
    left_join(rep_to_comp) %>%
  filter(Depth_Class == "E") %>% 
    st_as_sf(coords = c("Longitude", "Latitude"), crs = "EPSG:4326") %>%
    st_transform(crs = "EPSG:3174") %>%
  select(Rep_ID, AVG_W_DIST = avg.dist, AVG_UW_DIST = avg.un_dist, month) %>% # For some ungodly reason, QGIS HATES the avg.* name
  group_split(month)

st_write(split_dists[[1]], "analysis/QGIS_Work/May_Avg_Dists.gpkg" , append = FALSE)
## Deleting layer `May_Avg_Dists' using driver `GPKG'
## Writing layer `May_Avg_Dists' to data source `analysis/QGIS_Work/May_Avg_Dists.gpkg' using driver `GPKG'
## Writing 15 features with 4 fields and geometry type Point.
st_write(split_dists[[2]], "analysis/QGIS_Work/Sept_Avg_Dists.gpkg", append = FALSE)
## Deleting layer `Sept_Avg_Dists' using driver `GPKG'
## Writing layer `Sept_Avg_Dists' to data source `analysis/QGIS_Work/Sept_Avg_Dists.gpkg' using driver `GPKG'
## Writing 15 features with 4 fields and geometry type Point.

Now, we calculate other pairwise distances, including geographic distance and depth difference.

geo_dist_mat <- sample_sf %>%
  st_distance()

row.names(geo_dist_mat) <- sample_sf$Rep_ID
colnames(geo_dist_mat) <- sample_sf$Rep_ID

diag(geo_dist_mat) <- NA

geo_dist_mat[upper.tri(geo_dist_mat)] <- NA

geo_mat_df <- as.data.frame(geo_dist_mat)

tidy_geo_dist <- geo_mat_df %>%
    mutate(item1 = rownames(geo_mat_df)) %>%
    pivot_longer(!item1, names_to = "item2", values_to = "geo_dist") %>%
  filter(!is.na(geo_dist))
just_depth <- sample_sf %>%
  dplyr::select(Rep_ID, Deployment_Depth_m)

depth_dist_mat <- abs(outer(just_depth$Deployment_Depth_m, just_depth$Deployment_Depth_m, "-"))


row.names(depth_dist_mat) <- sample_sf$Rep_ID
colnames(depth_dist_mat) <- sample_sf$Rep_ID

diag(depth_dist_mat) <- NA
depth_dist_mat[upper.tri(depth_dist_mat)] <- NA

depth_mat_df <- as.data.frame(depth_dist_mat)

tidy_depth_dist <- depth_mat_df %>%
    mutate(item1 = rownames(depth_mat_df)) %>%
    pivot_longer(!item1, names_to = "item2", values_to = "depth_dist") %>%
  filter(!is.na(depth_dist))
for_join <- sample_sf %>%
  as.data.frame() %>%
  dplyr::select(Rep_ID, Comp_Group_Hier, month, -geometry)

stopifnot(all(tidy_geo_dist$item1 == tidy_depth_dist$item1))
stopifnot(all(tidy_geo_dist$item2 == tidy_depth_dist$item2))


tidy_dist_decay <- tidy_dist %>%
  full_join(tidy_geo_dist) %>%
  full_join(tidy_depth_dist) %>%
  filter(!is.na(geo_dist)) %>%
  left_join(for_join, by = c("item1" = "Rep_ID")) %>%
  left_join(for_join, by = c("item2" = "Rep_ID"))

Next, we make many different distance decays, based on either unweighted vs. weighted, and deep vs. shallow samples. Note that we make the linear models first and then use them to annotate the plots

geo_models <- tidy_dist_decay %>%
  filter(month.x == month.y,
         str_detect(Comp_Group_Hier.x, "Shallow"),
         str_detect(Comp_Group_Hier.y, "Shallow")) %>%
  mutate(geo_dist = as.numeric(geo_dist)) %>%
  nest_by(month.x) %>%
  mutate(lms = list(lm(1-dist ~ geo_dist, data = data))) %>%
  pull(lms)

labs <- map(geo_models, \(x){
  r <-summary(x)$adj.r.squared
  p <- summary(x)$coefficients[2,4]
  r.lab <- ifelse(r < 0.01,
        "R[adj]^2 < 0.01~~~",
        paste("R[adj]^2 ==", round(r,2), "~~~")
        )
  
  p.lab <- ifelse(p < 0.01,
                  "p < 0.01",
                  paste("p == ", round(p, 2))
                  )
  paste0(r.lab, p.lab)
  }
  )

dist_decay_weighted_surface <-
  tidy_dist_decay %>%
  filter(month.x == month.y,
         str_detect(Comp_Group_Hier.x, "Shallow"),
         str_detect(Comp_Group_Hier.y, "Shallow")) %>%
  mutate(geo_dist = as.numeric(geo_dist)) %>%
  ggplot(aes(x = geo_dist/1000, y = 1-dist, color = month.x)) + 
  geom_point(alpha = 0.2) + 
  scale_y_continuous(limits = c(0, 1)) + 
  labs(x = "Geographic Distance (km)", y = "Similarity (Weighted)", color = "Month") + 
  geom_smooth(method = "lm", se = FALSE) + 
  scale_color_manual(values = triglav[c(1,5)]) + 
  annotate(geom = "text",
           x = 0, 
           y = c(.33,1),
           color = triglav[c(1,5)],
           hjust = 0,
           label = labs,
           parse = TRUE)
geo_models <- tidy_dist_decay %>%
  filter(month.x == month.y,
         str_detect(Comp_Group_Hier.x, "Shallow"),
         str_detect(Comp_Group_Hier.y, "Shallow")) %>%
  mutate(geo_dist = as.numeric(geo_dist)) %>%
  nest_by(month.x) %>%
  mutate(lms = list(lm(1-un_dist ~ geo_dist, data = data))) %>%
  pull(lms)

labs <- map(geo_models, \(x){
  r <-summary(x)$adj.r.squared
  p <- summary(x)$coefficients[2,4]
  r.lab <- ifelse(r < 0.01,
        "R[adj]^2 < 0.01~~~",
        paste("R[adj]^2 ==", round(r,2), "~~~")
        )
  
  p.lab <- ifelse(p < 0.01,
                  "p < 0.01",
                  paste("p == ", round(p, 2))
                  )
  paste0(r.lab, p.lab)
  }
  )

dist_decay_unweighted_shallow <- 
  tidy_dist_decay %>%
  filter(month.x == month.y,
         str_detect(Comp_Group_Hier.x, "Shallow"),
         str_detect(Comp_Group_Hier.y, "Shallow")) %>%
  mutate(geo_dist = as.numeric(geo_dist)) %>%
  ggplot(aes(x = geo_dist/1000, y = 1-un_dist, color = month.x)) + 
  geom_point(alpha = 0.2) + 
  scale_y_continuous(limits = c(0, 1)) + 
  labs(x = "Geographic Distance (km)", y = "Similarity (Unweighted)", color = "Month") + 
  geom_smooth(method = "lm", se = FALSE) + 
  scale_color_manual(values = triglav[c(1,5)]) + 
  annotate(geom = "text",
           x = 0, 
           y = c(0.25,1),
           color = triglav[c(1,5)],
           hjust = 0,
           label = labs,
           parse = TRUE)
geo_models <- tidy_dist_decay %>%
  filter(month.x == month.y,
         Comp_Group_Hier.x == "Deep",
         Comp_Group_Hier.y == "Deep") %>%
  mutate(geo_dist = as.numeric(geo_dist)) %>%
  nest_by(month.x) %>%
  mutate(lms = list(lm(1-dist ~ geo_dist, data = data))) %>%
  pull(lms)

labs <- map(geo_models, \(x){
  r <-summary(x)$adj.r.squared
  p <- summary(x)$coefficients[2,4]
  r.lab <- ifelse(r < 0.01,
        "R[adj]^2 < 0.01~~~",
        paste("R[adj]^2 ==", round(r,2), "~~~")
        )
  
  p.lab <- ifelse(p < 0.01,
                  "p < 0.01",
                  paste("p == ", round(p, 2))
                  )
  paste0(r.lab, p.lab)
  }
  )

dist_decay_weighted_deep <- 
  tidy_dist_decay %>%
  filter(month.x == month.y,
         str_detect(Comp_Group_Hier.x, "Deep"),
         str_detect(Comp_Group_Hier.y, "Deep")) %>%
  mutate(geo_dist = as.numeric(geo_dist)) %>%
  ggplot(aes(x = geo_dist/1000, y = 1-dist, color = month.x)) + 
  geom_point(alpha = 0.2) + 
  scale_y_continuous(limits = c(0, 1)) + 
  labs(x = "Geographic Distance (km)", y = "Similarity (Weighted)", color = "Month") + 
  geom_smooth(method = "lm", se = FALSE) + 
  scale_color_manual(values = triglav[c(1,5)]) + 
  annotate(geom = "text",
           x = 0, 
           y = c(1,0.33),
           color = triglav[c(1,5)],
           hjust = 0,
           label = labs,
           parse = TRUE)
geo_models <- tidy_dist_decay %>%
  filter(month.x == month.y,
         str_detect(Comp_Group_Hier.x, "Deep"),
         str_detect(Comp_Group_Hier.y, "Deep")) %>%
  mutate(geo_dist = as.numeric(geo_dist)) %>%
  nest_by(month.x) %>%
  mutate(lms = list(lm(1-un_dist ~ geo_dist, data = data))) %>%
  pull(lms)

labs <- map(geo_models, \(x){
  r <-summary(x)$adj.r.squared
  p <- summary(x)$coefficients[2,4]
  r.lab <- ifelse(r < 0.01,
        "R[adj]^2 < 0.01~~~",
        paste("R[adj]^2 ==", round(r,2), "~~~")
        )
  
  p.lab <- ifelse(p < 0.01,
                  "p < 0.01",
                  paste("p == ", round(p, 2))
                  )
  paste0(r.lab, p.lab)
  }
  )

dist_decay_unweighted_deep <- 
  tidy_dist_decay %>%
  filter(month.x == month.y,
         str_detect(Comp_Group_Hier.x, "Deep"),
         str_detect(Comp_Group_Hier.y, "Deep")) %>%
  mutate(geo_dist = as.numeric(geo_dist)) %>%
  ggplot(aes(x = geo_dist/1000, y = 1-un_dist, color = month.x)) + 
  geom_point(alpha = 0.2) + 
  scale_y_continuous(limits = c(0, 1)) + 
  labs(x = "Geographic Distance (km)", y = "Similarity (Unweighted)", color = "Month") + 
  geom_smooth(method = "lm", se = FALSE) + 
  scale_color_manual(values = triglav[c(1,5)]) + 
  annotate(geom = "text",
           x = 0, 
           y = c(1,0.15),
           color = triglav[c(1,5)],
           hjust = 0,
           label = labs,
           parse = TRUE)
(dist_decay_unweighted_shallow + dist_decay_weighted_surface) / (dist_decay_unweighted_deep +dist_decay_weighted_deep) + 
  plot_layout(guides = "collect")

depth_models <- tidy_dist_decay %>%
  filter(month.x == month.y) %>%
  nest_by(month.x) %>%
  mutate(lms = list(lm(1-un_dist ~ depth_dist, data = data))) %>%
  pull(lms)

labs <- map(depth_models, \(x){
  r <-summary(x)$adj.r.squared
  p <- summary(x)$coefficients[2,4]
  r.lab <- ifelse(r < 0.01,
        "R[adj]^2 < 0.01~~~",
        paste("R[adj]^2 ==", round(r,2), "~~~")
        )
  
  p.lab <- ifelse(p < 0.01,
                  "p < 0.01",
                  paste("p == ", round(p, 2))
                  )
  paste0(r.lab, p.lab)
  }
  )

depth_decay_unweighted <- tidy_dist_decay %>%
  filter(month.x == month.y) %>%
  ggplot(aes(x = depth_dist, y = 1-un_dist, color = month.x)) + 
  geom_point(alpha = 0.2) + 
  scale_y_continuous(limits = c(0, 1)) + 
  labs(x = "Depth Difference (m)", y = "Similarity (Unweighted)", color = "Month") + 
  geom_smooth(method = "lm", se = FALSE) + 
  scale_color_manual(values = triglav[c(1,5)]) + 
  annotate(geom = "text",
           x = 0, 
           y = c(1,0),
           color = triglav[c(1,5)],
           hjust = 0,
           label = labs,
           parse = TRUE)


depth_decay_unweighted

Now we do the decay relationships for the main text, which use weighted distance and don’t separate by depths.

geo_models <- tidy_dist_decay %>%
  filter(month.x == month.y) %>%
  mutate(geo_dist = as.numeric(geo_dist)) %>%
  nest_by(month.x) %>%
  mutate(lms = list(lm(1-dist ~ geo_dist, data = data))) %>%
  pull(lms)

labs <- map(geo_models, \(x){
  r <-summary(x)$adj.r.squared
  p <- summary(x)$coefficients[2,4]
  r.lab <- ifelse(r < 0.01,
        "R[adj]^2 < 0.01~~~",
        paste("R[adj]^2 ==", round(r,2), "~~~")
        )
  
  p.lab <- ifelse(p < 0.01,
                  "p < 0.01",
                  paste("p == ", round(p, 2))
                  )
  paste0(r.lab, p.lab)
  }
  )

dist_decay <- tidy_dist_decay %>%
  filter(month.x == month.y) %>%
  mutate(geo_dist = as.numeric(geo_dist)) %>%
  ggplot(aes(x = geo_dist/1000, y = 1-dist, color = month.x)) + 
  geom_point(alpha = 0.2) + 
  scale_y_continuous(limits = c(0, 1)) + 
  labs(x = "Geographic Distance (km)", y = "Similarity", color = "Month") + 
  geom_smooth(method = "lm", se = FALSE) + 
  scale_color_manual(values = triglav[c(1,5)]) + 
  annotate(geom = "text",
           x = 0, 
           y = c(1,0.15),
           color = triglav[c(1,5)],
           hjust = 0,
           label = labs,
           parse = TRUE)

Now we do depth decay

depth_models <- tidy_dist_decay %>%
  filter(month.x == month.y) %>%
  nest_by(month.x) %>%
  mutate(lms = list(lm(1-dist ~ depth_dist, data = data))) %>%
  pull(lms)

labs <- map(depth_models, \(x){
  r <-summary(x)$adj.r.squared
  p <- summary(x)$coefficients[2,4]
  r.lab <- ifelse(r < 0.01,
        "R[adj]^2 < 0.01~~~",
        paste("R[adj]^2 ==", round(r,2), "~~~")
        )
  
  p.lab <- ifelse(p < 0.01,
                  "p < 0.01",
                  paste("p == ", round(p, 2))
                  )
  paste0(r.lab, p.lab)
  }
  )

depth_decay <- tidy_dist_decay %>%
  filter(month.x == month.y) %>%
  ggplot(aes(x = depth_dist, y = 1-dist, color = month.x)) + 
  geom_point(alpha = 0.2) + 
  scale_y_continuous(limits = c(0, 1)) + 
  labs(x = "Depth Difference (m)", y = "Similarity", color = "Month") + 
  geom_smooth(method = "lm", se = FALSE) + 
  scale_color_manual(values = triglav[c(1,5)]) + 
  annotate(geom = "text",
           x = 0, 
           y = c(1,0.15),
           color = triglav[c(1,5)],
           hjust = 0,
           label = labs,
           parse = TRUE)
combined <- ggpubr::ggarrange(dist_decay, depth_decay, 
                              nrow = 2,
                  common.legend = TRUE,
                  legend = "bottom")

combined

Next, we run variance partitioning, including testing for significance.

full_records <- full_abs_physeq %>%
  sample_data() %>%
  data.frame() %>%
  dplyr::select(par, temperature, TN, TP, Na, SO4, DOC, chl_a, Rep_ID) %>%
  drop_na() %>%
  pull(Rep_ID)

no_na_physeq <- full_abs_physeq %>%
  microViz::ps_filter(Rep_ID %in%full_records)

metadata <- no_na_physeq %>%
  sample_data() %>%
  data.frame()

month <- metadata["month"]
depth <- metadata["Deployment_Depth_m"]
geo <- dplyr::select(metadata, Longitude, Latitude)
env <- dplyr::select(metadata, NH4:chl_a,temperature,good_oxygen,par) %>% scale()

missing <- rownames(abs_unifrac_dist_0.5)[!rownames(abs_unifrac_dist_0.5) %in% rownames(month)]

missing
## [1] "May_41_M"
abs_wun <- as.dist(abs_unifrac_dist_0.5[setdiff(rownames(abs_unifrac_dist_0.5),missing),setdiff(rownames(abs_unifrac_dist_0.5),missing)])
abs_un <- as.dist(unweighted_unifrac_dist[setdiff(rownames(unweighted_unifrac_dist),missing),setdiff(rownames(unweighted_unifrac_dist),missing)])

metadata <- no_na_physeq %>%
  sample_data() %>%
  data.frame()



env_variables <- paste(colnames(env), collapse = "+")

env_form <- as.formula(paste0("abs_wun ~ ", env_variables))

env_d_form <- as.formula(paste0("abs_wun ~ ", env_variables, "+ Condition(Deployment_Depth_m)"))

env_m_form <- as.formula(paste0("abs_wun ~ ", env_variables, "+ Condition(month)"))

env_g_form <- as.formula(paste0("abs_wun ~ ", env_variables, "+ Condition(Latitude + Longitude)"))

m_form <- as.formula(paste0("abs_wun ~ ", "month"))

d_form <- as.formula(paste0("abs_wun ~ ", "Deployment_Depth_m"))

g_form <- as.formula(paste0("abs_wun ~ ", "Latitude + Longitude"))

formulas <- list(Env = env_form, 
                 Env_x_Depth = env_d_form, 
                 Env_x_Month = env_m_form, 
                 Env_x_Geo = env_g_form, 
                 Month = m_form,
                 Depth = d_form,
                 Geo = g_form)

varpart_results <- vegan::varpart(Y = abs_wun, month, depth, geo, env)

plot(varpart_results, Xnames = c("Month","Depth","Geo. Distance","Env."))

# Test the significance
map(formulas, \(form){
  dbrda <- vegan::dbrda(form, data = as.data.frame(cbind(env, geo, depth, month)))
  cca_res <- vegan::anova.cca(dbrda)
  cca_res$`Pr(>F)`[1]
})
## $Env
## [1] 0.001
## 
## $Env_x_Depth
## [1] 0.001
## 
## $Env_x_Month
## [1] 0.001
## 
## $Env_x_Geo
## [1] 0.001
## 
## $Month
## [1] 0.001
## 
## $Depth
## [1] 0.001
## 
## $Geo
## [1] 0.061

Then, we remake this figure, to be simpler and clearer:

vars <- c("Env.", 
          "Env. + Depth", 
          "Env. + Month", 
          "Depth", 
          "Month", 
          "Env. + Location",
          "Location",
          "Residuals")

var_df <- data.frame(Label = factor(vars, levels = rev(vars)),
           Variance = c(.24,
                        .22,
                        .15,
                        .03,
                        .03,
                        .02,
                        .01,
                        .32),
           Sig = factor(c(1,1,1,1,1,1,0,1)))

var_plot <- ggplot(var_df, 
       aes(y = Label,
           x = Variance,
           color = Sig)) + 
  geom_point() +
  geom_segment(aes(xend = 0, yend= Label)) + 
  scale_x_continuous(expand = expansion(mult = c(0,.15))) + 
  geom_text(aes(label = Variance,
                x = Variance + .01),
            hjust = 0) + 
  scale_color_manual(values = c("grey70","black"), guide = "none") + 
  theme(axis.title.y = element_blank())

var_plot

Then, we do variance partitioning with unweighted Unifrac.

full_records <- full_abs_physeq %>%
  sample_data() %>%
  data.frame() %>%
  dplyr::select(par, temperature, TN, TP, Na, SO4, DOC, chl_a, Rep_ID) %>%
  drop_na() %>%
  pull(Rep_ID)

env_form <- as.formula(paste0("abs_wun ~ ", env_variables))

env_d_form <- as.formula(paste0("abs_wun ~ ", env_variables, "+ Condition(Deployment_Depth_m)"))

env_m_form <- as.formula(paste0("abs_wun ~ ", env_variables, "+ Condition(month)"))

env_g_form <- as.formula(paste0("abs_wun ~ ", env_variables, "+ Condition(Latitude + Longitude)"))

m_form <- as.formula(paste0("abs_wun ~ ", "month"))

d_form <- as.formula(paste0("abs_wun ~ ", "Deployment_Depth_m"))

formulas <- list(Env = env_form, 
                 Env_x_Depth = env_d_form, 
                 Env_x_Month = env_m_form, 
                 Env_x_Geo = env_g_form, 
                 Month = m_form,
                 Depth = d_form)


varpart_results <- vegan::varpart(Y = abs_un, month, depth, geo, env)

plot(varpart_results, Xnames = c("Month","Depth","Geo. Distance","Env."))

I love this plot so I’m keeping it in here, but we don’t use it in the paper.

mat_df_wun <- as.data.frame(abs_unifrac_dist_0.5)

mat_df_un <- as.data.frame(unweighted_unifrac_dist)

tidy_dist_wun <- mat_df_wun %>%
    mutate(item1 = rownames(mat_df_wun)) %>%
    pivot_longer(!item1, names_to = "item2", values_to = "wunifrac") %>%
  filter(!is.na(wunifrac))

tidy_dist_un <- mat_df_un %>%
    mutate(item1 = rownames(mat_df_un)) %>%
    pivot_longer(!item1, names_to = "item2", values_to = "unifrac") %>%
  filter(!is.na(unifrac))

both_dists <- full_join(tidy_dist_wun, tidy_dist_un)

metadata <- no_na_physeq %>%
  sample_data() %>%
  data.frame()

env <- dplyr::select(metadata, NH4:chl_a,temperature,good_oxygen,par) %>% scale()

tidy_euc <- 
  env %>%
  dist() %>% 
  as.matrix() %>% 
  data.frame() %>% 
  mutate(item1 = rownames(.)) %>% 
  mutate(across(everything(), \(x)ifelse(x == 0, NA, x))) %>%
  pivot_longer(!item1, names_to = "item2", values_to = "euc_dist") %>%
  filter(!is.na(euc_dist))

wun_and_euc <- full_join(both_dists, tidy_euc) %>%
  full_join(tidy_geo_dist) %>%
  filter(!is.na(geo_dist)) %>%
  transmute(across(where(is.numeric), scale))



wun_and_euc %>%
  pivot_longer(c(wunifrac,unifrac), names_to = "Community_Metric", values_to = "Community_Distance") %>%
  pivot_longer(c(euc_dist, geo_dist), names_to = "Independent_Metric",values_to = "Independent_Distance") %>%
  ggplot(aes(x = Independent_Distance, y = Community_Distance)) + 
  geom_point(alpha = 0.1) + 
  geom_smooth(method = "lm", se = FALSE) +
  facet_grid(Independent_Metric ~ Community_Metric)

wun_and_euc %>%
  select(-unifrac) %>%
  pivot_longer(euc_dist:geo_dist, names_to = "Metric", values_to = "Distance") %>%
  ggplot(aes(x = Distance, y = wunifrac, color = Metric)) + 
  geom_point(alpha = 0.1) + 
  geom_smooth(method = "lm", se = FALSE, linewidth = 2,
              linetype = 2) + 
  labs(x = "Scaled Geographic or Environmental Distance", y = "Weighted Unifrac Distance",
       color = "Distance") + 
  scale_color_manual(labels = c("Environmental","Geographical"),
                     values = triglav[c(1,4)])

Session Info

sessioninfo::session_info()
## ─ Session info ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.3.3 (2024-02-29)
##  os       Rocky Linux 9.5 (Blue Onyx)
##  system   x86_64, linux-gnu
##  ui       X11
##  language (EN)
##  collate  en_US.UTF-8
##  ctype    en_US.UTF-8
##  tz       America/New_York
##  date     2025-06-27
##  pandoc   3.1.1 @ /usr/lib/rstudio-server/bin/quarto/bin/tools/ (via rmarkdown)
## 
## ─ Packages ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
##  ! package           * version     date (UTC) lib source
##  P abind               1.4-5       2016-07-21 [?] CRAN (R 4.3.2)
##  P ade4                1.7-22      2023-02-06 [?] CRAN (R 4.3.2)
##  P ape                 5.7-1       2023-03-13 [?] CRAN (R 4.3.2)
##  P backports           1.4.1       2021-12-13 [?] CRAN (R 4.3.2)
##  P base64enc           0.1-3       2015-07-28 [?] CRAN (R 4.3.2)
##  P beeswarm            0.4.0       2021-06-01 [?] CRAN (R 4.3.2)
##  P Biobase             2.62.0      2023-10-24 [?] Bioconductor
##  P BiocGenerics        0.48.1      2023-11-01 [?] Bioconductor
##  P BiocManager         1.30.22     2023-08-08 [?] CRAN (R 4.3.2)
##  P biomformat          1.30.0      2023-10-24 [?] Bioconductor
##  P Biostrings          2.70.1      2023-10-25 [?] Bioconductor
##  P bitops              1.0-7       2021-04-24 [?] CRAN (R 4.3.2)
##  P broom               1.0.5       2023-06-09 [?] CRAN (R 4.3.2)
##  P bslib               0.5.1       2023-08-11 [?] CRAN (R 4.3.2)
##  P cachem              1.0.8       2023-05-01 [?] CRAN (R 4.3.2)
##  P car                 3.1-2       2023-03-30 [?] CRAN (R 4.3.2)
##  P carData             3.0-5       2022-01-06 [?] CRAN (R 4.3.2)
##  P class               7.3-22      2023-05-03 [?] CRAN (R 4.3.3)
##  P classInt            0.4-10      2023-09-05 [?] CRAN (R 4.3.2)
##  P cli                 3.6.1       2023-03-23 [?] CRAN (R 4.3.2)
##  P cluster             2.1.4       2022-08-22 [?] CRAN (R 4.3.2)
##  P codetools           0.2-19      2023-02-01 [?] CRAN (R 4.3.3)
##  P colorspace          2.1-0       2023-01-23 [?] CRAN (R 4.3.2)
##  P cols4all            0.7-1       2024-03-12 [?] CRAN (R 4.3.3)
##  P cowplot             1.1.3       2024-01-22 [?] CRAN (R 4.3.2)
##  P crayon              1.5.2       2022-09-29 [?] CRAN (R 4.3.2)
##  P crosstalk           1.2.1       2023-11-23 [?] CRAN (R 4.3.2)
##  P cytolib             2.14.1      2024-01-18 [?] Bioconduc~
##  P data.table          1.15.2      2024-02-29 [?] CRAN (R 4.3.2)
##  P DBI                 1.2.2       2024-02-16 [?] CRAN (R 4.3.2)
##  P dichromat           2.0-0.1     2022-05-02 [?] CRAN (R 4.3.2)
##  P digest              0.6.33      2023-07-07 [?] CRAN (R 4.3.2)
##  P dplyr             * 1.1.3       2023-09-03 [?] CRAN (R 4.3.2)
##  P e1071               1.7-14      2023-12-06 [?] CRAN (R 4.3.2)
##  P ellipsis            0.3.2       2021-04-29 [?] CRAN (R 4.3.2)
##  P evaluate            0.23        2023-11-01 [?] CRAN (R 4.3.2)
##  P fansi               1.0.5       2023-10-08 [?] CRAN (R 4.3.2)
##  P farver              2.1.1       2022-07-06 [?] CRAN (R 4.3.2)
##  P fastmap             1.1.1       2023-02-24 [?] CRAN (R 4.3.2)
##  P flowCore            2.14.2      2024-03-18 [?] Bioconduc~
##  P forcats           * 1.0.0       2023-01-29 [?] CRAN (R 4.3.2)
##  P foreach             1.5.2       2022-02-02 [?] CRAN (R 4.3.2)
##  P fun.gus           * 0.3.1       2025-06-26 [?] Github (MarschmiLab/fun.gus@7daa3fa)
##  P furrr               0.3.1       2022-08-15 [?] CRAN (R 4.3.2)
##  P future              1.33.1      2023-12-22 [?] CRAN (R 4.3.3)
##  P generics            0.1.3       2022-07-05 [?] CRAN (R 4.3.2)
##  P GenomeInfoDb        1.38.0      2023-10-24 [?] Bioconductor
##  P GenomeInfoDbData    1.2.11      2023-11-07 [?] Bioconductor
##  P ggbeeswarm          0.7.2       2023-04-29 [?] CRAN (R 4.3.2)
##  P ggplot2           * 3.5.0       2024-02-23 [?] CRAN (R 4.3.2)
##  P ggpubr              0.6.0       2023-02-10 [?] CRAN (R 4.3.2)
##  P ggsignif            0.6.4       2022-10-13 [?] CRAN (R 4.3.2)
##  P globals             0.16.2      2022-11-21 [?] CRAN (R 4.3.3)
##  P glue                1.6.2       2022-02-24 [?] CRAN (R 4.3.2)
##  P gridExtra           2.3         2017-09-09 [?] CRAN (R 4.3.2)
##  P gtable              0.3.4       2023-08-21 [?] CRAN (R 4.3.2)
##  P highr               0.10        2022-12-22 [?] CRAN (R 4.3.2)
##  P hms                 1.1.3       2023-03-21 [?] CRAN (R 4.3.2)
##  P htmltools           0.5.7       2023-11-03 [?] CRAN (R 4.3.2)
##  P htmlwidgets         1.6.2       2023-03-17 [?] CRAN (R 4.3.2)
##  P httpuv              1.6.12      2023-10-23 [?] CRAN (R 4.3.2)
##  P igraph              1.5.1       2023-08-10 [?] CRAN (R 4.3.2)
##  P IRanges             2.36.0      2023-10-24 [?] Bioconductor
##  P isoband             0.2.7       2022-12-20 [?] CRAN (R 4.3.2)
##  P iterators           1.0.14      2022-02-05 [?] CRAN (R 4.3.2)
##  P jquerylib           0.1.4       2021-04-26 [?] CRAN (R 4.3.2)
##  P jsonlite            1.8.7       2023-06-29 [?] CRAN (R 4.3.2)
##  P KernSmooth          2.23-22     2023-07-10 [?] CRAN (R 4.3.3)
##  P knitr               1.45        2023-10-30 [?] CRAN (R 4.3.2)
##  P labeling            0.4.3       2023-08-29 [?] CRAN (R 4.3.2)
##  P later               1.3.1       2023-05-02 [?] CRAN (R 4.3.2)
##  P lattice             0.21-9      2023-10-01 [?] CRAN (R 4.3.2)
##  P leafem              0.2.3       2023-09-17 [?] CRAN (R 4.3.2)
##  P leaflegend          1.2.0       2024-01-10 [?] CRAN (R 4.3.3)
##  P leaflet             2.2.1       2023-11-13 [?] CRAN (R 4.3.2)
##  P leaflet.providers   2.0.0       2023-10-17 [?] CRAN (R 4.3.2)
##  P leafsync            0.1.0       2019-03-05 [?] CRAN (R 4.3.2)
##  P lifecycle           1.0.3       2022-10-07 [?] CRAN (R 4.3.2)
##  P listenv             0.9.1       2024-01-29 [?] CRAN (R 4.3.2)
##  P lubridate         * 1.9.3       2023-09-27 [?] CRAN (R 4.3.2)
##  P lwgeom              0.2-14      2024-02-21 [?] CRAN (R 4.3.2)
##  P magrittr            2.0.3       2022-03-30 [?] CRAN (R 4.3.2)
##  P MASS                7.3-60      2023-05-04 [?] CRAN (R 4.3.2)
##  P Matrix              1.6-1.1     2023-09-18 [?] CRAN (R 4.3.2)
##  P matrixStats         1.2.0       2023-12-11 [?] CRAN (R 4.3.2)
##  P MBA                 0.1-0       2022-11-29 [?] CRAN (R 4.3.2)
##  P mgcv                1.9-0       2023-07-11 [?] CRAN (R 4.3.2)
##  P microViz            0.12.1      2024-03-05 [?] https://d~
##  P mime                0.12        2021-09-28 [?] CRAN (R 4.3.2)
##  P multtest            2.58.0      2023-10-24 [?] Bioconductor
##  P munsell             0.5.0       2018-06-12 [?] CRAN (R 4.3.2)
##  P NatParksPalettes  * 0.2.0       2022-10-09 [?] CRAN (R 4.3.2)
##  P nlme                3.1-163     2023-08-09 [?] CRAN (R 4.3.2)
##  P pacman              0.5.1       2019-03-11 [?] CRAN (R 4.3.2)
##  P parallelly          1.37.1      2024-02-29 [?] CRAN (R 4.3.3)
##  P patchwork         * 1.2.0.9000  2025-06-26 [?] Github (thomasp85/patchwork@d943757)
##  P permute             0.9-7       2022-01-27 [?] CRAN (R 4.3.2)
##  P phyloseq          * 1.46.0      2023-10-24 [?] Bioconductor
##  P pillar              1.9.0       2023-03-22 [?] CRAN (R 4.3.2)
##  P pkgconfig           2.0.3       2019-09-22 [?] CRAN (R 4.3.2)
##  P plyr                1.8.9       2023-10-02 [?] CRAN (R 4.3.2)
##  P png                 0.1-8       2022-11-29 [?] CRAN (R 4.3.2)
##  P promises            1.2.1       2023-08-10 [?] CRAN (R 4.3.2)
##  P proxy               0.4-27      2022-06-09 [?] CRAN (R 4.3.2)
##  P purrr             * 1.0.2       2023-08-10 [?] CRAN (R 4.3.2)
##  P R6                  2.5.1       2021-08-19 [?] CRAN (R 4.3.2)
##  P raster              3.6-26      2023-10-14 [?] CRAN (R 4.3.2)
##  P RColorBrewer        1.1-3       2022-04-03 [?] CRAN (R 4.3.2)
##  P Rcpp                1.0.11      2023-07-06 [?] CRAN (R 4.3.2)
##  P RCurl               1.98-1.13   2023-11-02 [?] CRAN (R 4.3.2)
##  P readr             * 2.1.5       2024-01-10 [?] CRAN (R 4.3.2)
##    renv                1.0.5       2024-02-29 [1] CRAN (R 4.3.2)
##  P reshape2            1.4.4       2020-04-09 [?] CRAN (R 4.3.2)
##  P rhdf5               2.46.1      2023-11-29 [?] Bioconduc~
##  P rhdf5filters        1.14.1      2023-11-06 [?] Bioconductor
##  P Rhdf5lib            1.24.2      2024-02-07 [?] Bioconduc~
##  P rlang               1.1.2       2023-11-04 [?] CRAN (R 4.3.2)
##  P rmarkdown           2.25        2023-09-18 [?] CRAN (R 4.3.2)
##  P RProtoBufLib        2.14.1      2024-03-18 [?] Bioconduc~
##  P rstatix             0.7.2       2023-02-01 [?] CRAN (R 4.3.2)
##  P rstudioapi          0.15.0      2023-07-07 [?] CRAN (R 4.3.2)
##  P s2                  1.1.6       2023-12-19 [?] CRAN (R 4.3.2)
##  P S4Vectors           0.40.1      2023-10-26 [?] Bioconductor
##  P sass                0.4.7       2023-07-15 [?] CRAN (R 4.3.2)
##  P scales              1.3.0       2023-11-28 [?] CRAN (R 4.3.2)
##  P sessioninfo         1.2.2       2021-12-06 [?] CRAN (R 4.3.2)
##  P sf                * 1.0-15      2023-12-18 [?] CRAN (R 4.3.2)
##  P shiny               1.7.5.1     2023-10-14 [?] CRAN (R 4.3.2)
##  P sp                  2.1-3       2024-01-30 [?] CRAN (R 4.3.3)
##  P spacesXYZ           1.3-0       2024-01-23 [?] CRAN (R 4.3.3)
##  P stars               0.6-4       2023-09-11 [?] CRAN (R 4.3.2)
##  P stringi             1.7.12      2023-01-11 [?] CRAN (R 4.3.2)
##  P stringr           * 1.5.0       2022-12-02 [?] CRAN (R 4.3.2)
##  P survival            3.5-8       2024-02-14 [?] CRAN (R 4.3.3)
##  P terra             * 1.7-71      2024-01-31 [?] CRAN (R 4.3.2)
##  P tibble            * 3.2.1       2023-03-20 [?] CRAN (R 4.3.2)
##  P tidyr             * 1.3.1       2024-01-24 [?] CRAN (R 4.3.2)
##  P tidyselect          1.2.0       2022-10-10 [?] CRAN (R 4.3.2)
##  P tidyverse         * 2.0.0       2023-02-22 [?] CRAN (R 4.3.2)
##  P timechange          0.3.0       2024-01-18 [?] CRAN (R 4.3.2)
##  P tmap              * 3.99.9000   2025-06-26 [?] Github (r-tmap/tmap@bd12c21)
##  P tmaptools         * 3.1-1       2025-06-26 [?] Github (mtennekes/tmaptools@0c8b0b1)
##  P tzdb                0.4.0       2023-05-12 [?] CRAN (R 4.3.2)
##  P units             * 0.8-5       2023-11-28 [?] CRAN (R 4.3.2)
##  P utf8                1.2.4       2023-10-22 [?] CRAN (R 4.3.2)
##  P vctrs               0.6.4       2023-10-12 [?] CRAN (R 4.3.2)
##  P vegan               2.6-4       2022-10-11 [?] CRAN (R 4.3.2)
##  P vipor               0.4.7       2023-12-18 [?] CRAN (R 4.3.2)
##  P viridisLite         0.4.2       2023-05-02 [?] CRAN (R 4.3.2)
##  P widgetframe         0.3.1       2017-12-20 [?] CRAN (R 4.3.2)
##  P withr               2.5.2       2023-10-30 [?] CRAN (R 4.3.2)
##  P wk                  0.9.1       2023-11-29 [?] CRAN (R 4.3.2)
##  P xfun                0.52        2025-04-02 [?] CRAN (R 4.3.3)
##  P XML                 3.99-0.16.1 2024-01-22 [?] CRAN (R 4.3.2)
##  P xtable              1.8-4       2019-04-21 [?] CRAN (R 4.3.2)
##  P XVector             0.42.0      2023-10-24 [?] Bioconductor
##  P yaml                2.3.7       2023-01-23 [?] CRAN (R 4.3.2)
##  P zlibbioc            1.48.0      2023-10-24 [?] Bioconductor
## 
##  [1] /local/workdir/arp277/Pendleton_2025_Ontario_Publication_Repo/renv/library/R-4.3/x86_64-pc-linux-gnu
##  [2] /home/arp277/.cache/R/renv/sandbox/R-4.3/x86_64-pc-linux-gnu/fd835031
## 
##  P ── Loaded and on-disk path mismatch.
## 
## ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────